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PREFACE 


This  report  was  prepared  by  Dr.  Giles  M.  Marion,  Research  Physical  Scientist,  Geochemical 
Sciences  Branch,  Research  Division,  U.S.  Army  Cold  Regions  Research  and  Engineering  Laboratory. 
Funding  was  provided  by  DA  Project  4A161 102AT24,  Research  in  Snow,  Ice  and  Frozen  Ground', 
Task  SS;  Work  Unit  020,  Transport  of  Chemical  Species  in  Soils. 

The  author  thanks  Dean  Pidgeon  for  his  assistance  in  developing  the  MS-DOS  version  of  the 
lONPAIR  program.  Technical  review  was  provided  by  Dr.  C.M.  Reynolds,  D.C.  Leggett  and  Dr.  I.K. 
Iskandar,  all  of  CRREL. 

Anyone  wishing  either  a  Macintosh  or  MS-DOS  version  of  this  program  on  disk  should  contact  Dr 
Marion  by  writing  to  him  at  CECRL-RC,  72  Lyme  Road,  Hanover  New  Hampshire  03755-1290,  or 
by  phoning  him  at  603-646-4676. 

The  contents  of  this  report  are  not  to  be  used  for  advertising  or  promotional  purposes.  Citation  of 
brand  names  does  not  constitute  an  official  endorsement  or  approval  of  the  use  of  such  commercial 
products. 


lONPAIR:  A  Chemical  Speciation  Program 
for  Calcareous  and  Gypsiferous  Soil  Solutions 

GILES  M.  MARION 


INTRODUCTION 

Calcareous  and  gypsiferous  soils  are  common  along 
many  river  floodplains  in  the  subarctic  and  arctic  regions 
of  Alaska  (Marion  et  al. ,  in  press  b)  as  well  as  in  more  arid 
regions  such  as  deserts  and  grasslands  (Schlesinger 
1985;  Marion  etal.  1990;  Marion  etal.,  in  press  a;  Reddy 
et  al.  1990;  Suarez  et  al.  1990).  The  lONPAIR  program 
was  designed  to  speciate  the  chemical  composition  of 
calcareous  and  gypsiferous  aqueous  solutions.  There  are 
a  number  of  commercially  available  chemical  speciation 
programs.  However,  some  programs  such  as  WATEQ4F 
(Ball  and  Nordstrom  1987)  do  not  allow  Pco2  0^'^ 
partial  pressure  of  CO2  [g])  as  an  input,  while  other 
programs  such  as  GEOCHEM  handle  PCO2  ^  incon¬ 
venient  fashion  requiring  a  trial-and-enror  approach  to 
solution  (Sposito  and  Mattigod  1979,  Suarez  1990). 
Since 
studie 

specify  Pco^  would  be  highly  useful,  especially  in 
laboratory  studies  of  calcareous  soils.  In  contrast  to  most 
existing  chemical  speciation  programs,  lONPAIR  was 
developed  as  a  "stand-alone”  program,  which  means  that 
one  does  not  need  a  programming  software  package  such 
as  BASIC  or  FORTRAN  in  order  to  run  the  program. 
Also,  in  contrast  to  most  other  programs,  lONPAIR  was 
developed  in  both  Macintosh  and  MS-DOS  (IBM-com¬ 
patible)  versions. 

The  objective  of  this  report  is  to  document  the 
lONPAIR  program.  Important  features  of  the  lONPAIR 
program  are  variable  inputs,  accuracy,  portability  and 
easy  utilization. 


PcOj  is  ^  easily  controlled  variable  in  laboratory 
s,  a  chemical  speciation  program  that  allows  one  to 


PROGRAM  DESCRIPTION 

Chemical  speciation  programs  partition  the  total  (mea¬ 
sured)  concentration  of  aqueous-phase  constitaents  into 


the  concentrations  of  individual  chemical  species  (the 
unknowns).  For  every  unknown  aqueous-phase  con¬ 
stituent,  an  independent  equation  is  needed.  In  the 
lONPAIR  program,  there  are  either  17  or  19  unknowns 
(or  equations). 

The  program  calculates  aqueous-phase  concentra¬ 
tions  by  successive  approximations  using  the  Newton- 
Raphson  algorithm  to  solve  a  set  of  non-linearequations 
(Kom  and  Kom  1968).  The  requisite  set  of  equations 
consists  of  mass  balance  equations  for  the  major  con¬ 
stituents  (Ca,  Mg,  K,  Na,  SO4  and  alkalinity),  which  are 
defined  in  terms  of  concentrations,  for  example 

Car  =  [Ca^]  4-  [CaHCO.^-^] 

-t- [CaCO30] -I- [CaS04‘^)  (1) 

and  stability  relations  forthe  ion-pairs  ai.d  carbonic  acid 
(Table  1),  which  are  defined  in  terp  .  of  activities,  for 
example 

(Ca-^-^)  (SO4-  -)/(CaS04'J)  =  K  (2) 

where  brackets  refer  to  concentrations  and  parentheses 
refer  to  activities. 

In  order  to  so’vp  the  non-linear  equations,  one  must 
estimate  the  a  iivity  coefficients  (y,)  that  are  related  to 
concentrations  (nij)  and  activities  (a;)  by 

(3) 

where  the  subscript  i  refers  to  the  /’*'  constituent  in 
solution. 

The  activity  coefficients  of  monovalent  and  divalent 
constituents  are  estimated  with  the  Davies  equation 

log(7,)  =  -A;,M[V7/(i.0+  V7)]  -  0.3/}  (4) 


Table  1.  The  pK  values  at  25°C  for  tiic  solution 
phase  complexes  between  cation  (C)  and  anion 
(A),  where  pK  =  -log  (C  x  \/CA  complex).  The 
equilibrium  constants  wc.  e  calculated  from  standard 
free  energy  data  (I  indsay  1979).  In  addition  to  the 
above  constants,  the  program  also  uses  the  Henry’s 
Law  constant  for  CO,  equilibrium  between  the  aque¬ 
ous  and  solution  phases  (pK  =  1.464). 


Anion 

Ca 

w.e 

CiUmn 

K 

N(i 

It 

SOj 

2..t09 

2.229 

o.x.so 

0.696 

HCO^ 

1.129 

1 .070 

-0.249* 

-0.249 

6.3W 

co' 

3.1.S2 

3.240 

1.268* 

1.268 

10.330 

•  .Assumed  lo  be  the  same  as  the  Na  analogues. 


where  A  is  aconstant  (=  0.5 1 1 6  at  25°C),  r,  is  the  valence 
of  the  /’*'  ion,  and  /  is  the  ionic  strength 

/  =  0.5  I  Wi (5) 

(Davies  1962).  Neutrally  charged  con.stituents  (e.g.. 
CaS04'^)  are  assumed  to  have  activity  coefficients  of 
1.00.  For  the  justification  of  this  assumption,  see  the 
discu.s.sion  in  the  Limitations  section. 

In  eq  1,  Cay  is  the  measured  total  Ca  concentration 
and  the  constituents  on  the  right  are  the  unknown  con¬ 
centrations.  The  activities  in  eq  2  must  al.so  be  e.xpressed 
rn  terms  of  concentrations  by  using  the  relation  in  eq  3. 
for  example 

ICa^^llSOa-'l/lCaSOV’l  =  /^/yca-YsOj  =  K’ .  (6) 

Initially,  the  concentrations  (m\),  ionic  strength  (/) 
and  activity  coefficients  (y,)  are  unknown.  Rough  esti¬ 
mates  of  the  unknown  concentrations  are  made  that 
allow  one  to  estimate  the  ionic  strength  (eq  5)  and  the 
activity  coefficients  (eq  4).  Then,  the  initial  estimates  of 
the  unknown  concentrations  are  refined  using  the  New- 
ton-Raphson  algorithm.  The  new  estimates  of  the  con¬ 
centrations  are  then  used  to  re-estimate  the  ionic  strength 
and  activity  coefficients,  which  in  turn  are  used  to  re- 
estimate  the  concentrations  with  the  Newton  -Raphson 
algorithm.  This  interactive  process  continues  until  con¬ 
secutive  estimates  of  every  unknow  n  concentration  agree 
to  within  ±0.1  . 

The  computer  program  ( Appendix  A)  has  three  input 
options  requiring  aspecification  of  1 )  pH  and  alkalinity. 
2)  P('Q,  atxi  alkalinity,  or  3)  Pc  (),and  pH.  Given  any  one 
of  these  options,  plus  measurements  of  total  dissolved 
Ca.  Mg.  K.  Na.  Cl  and  SO4.  the  program  calculates  the 
true  ionic  concentrations  in  the  solution  phase.  Ion-pairs 


are  considered  betw  een  the  cations.  CYC’',  Mg’^,  K  ’  and 
Na"’'.  and  the  anions.  HCO,”.  CO,“  and  S04"  ■^.  Ion- 
pairs  with  Cl '  are  ignored  because  they  only  form  weak 
complexes  (Lindsay  1979). 

Input  to  the  program  (Appendix  B)  is  via  keyboard 
from  computer  screen  queries.  First,  a  “Title”  for  the  data 
set  is  requested:  this  can  be  any  alphanumeric  statement 
up  to  80  characters.  Next,  you  must  choose  the  input 
option  (.see  previous  paragraph).  Then,  total  chemical 
concentrations  of  each  constituent  are  requested.  Except 
for  pH  (no  units),  alkalinity  (mol^^^  L"';  where  mol^  = 
equivalents),  and  Pco,  (turn),  the  input  concentration 
units  are  mol  L"' .  The  latter  are  most  conveniently  input 
in ■■e-fonnat”(i.e.. 0.0107  mol  L"'  =  1.07e-2,  Appendix 
B). 

Output  from  the  program  (Appendix  B)  includes  a 
charge  balance,  the  concentrations  and  activities  of  the 
individual  chemical  species,  the  ionic  strength,  and 
plAP*(CaC03)  and  plAPICaSOa)  where 

plAP(CaC03)  =  -log|(Ca-’-’')(CO.r-))  .  (7) 

Calculated  pi  APfCaCOx)  and  pi  AP(CaS04)  can  be  com¬ 
pared  to  the  solubility  products  forcalcite  (/?Arsp  =  8.48) 
and  gypsum  {pK^^  =  4.64)  respectively. 

The  computer  program  is  written  in  TrueBASIC  for 
both  Macintosh  and  MS-DOS  computers;  there  are  two 
versions  for  each  operating  system.'  One  version  re¬ 
quires  the  TrueBASIC  software  package  to  open  and 
run;  given  the  latter  soft  ware,  this  version  can  be  changed 
by  the  user  (e.g.,  input  and  output  fonnats.  stability 
constants,  temperature  adjustments,  etc.).  The  .second 
version  is  a  "stand-alone”  (or  executable  file)  program 
that  does  not  require  the  TrueBASIC  software  package; 
however,  this  version  can  not  be  user-modified. 

This  program  was  initially  written  in  Fortran  in  1 970 
and  has  been  tested  w  ith  many  data  sets  over  the  years 
(Marion  and  Babcock  1977;  Schlesinger  1985;  Marion 
et  al.  1990;  Marion  et  al..  in  press  a.b).  However,  the 
author  accepts  no  responsibility  for  the  accuracy  of  these 
calculations  and  urges  users  to  verify  accuracy  for  them¬ 
selves.  Moreover,  there  is  no  assurance  that  the  New  ton  - 
Raphson  algorithm  will  always  converge.  If  problems 
develop  with  non-converging  data  sets,  then  program 
lines  1590- 1860  may  need  to  be  modified.  For  problems 
w  ith  the  program,  you  can  phone  the  author. 

*  Ion  activity  product. 

If  you  w  ant  a  copy  of  the  KDNP.AIR  program  on  disk,  please 
contact  Dr.  Giles  Marion.  L'.S.ACECRl  -RC,  72  I  ynie  Road. 
Hanover  New  Hampshire  (I37.S.S  and  specify  it  you  need  a 
Macintosh  version  01  an  M.S-DOS  version,  or  call  Dr.  Marion 
al  Wl.s  fvf6-4()76 


The  current  Macintosh  version  has  run  on  the  Macin¬ 
tosh  Plus,  SE,  Classic  and  II  computers  with  both 
ImageWriter  and  LaserWriter  printers.  The  MS-DOS 
version  has  also  been  tested  on  a  wide  range  of  MS-DOS 
computers. 

LIMITATIONS 

This  program  requires  a  specification  of  all  constitu¬ 
ents;  you  can  not  enter  a  "zero"  concentration.  The  latter 
will  cause  program  failure.  If  a  constituent  is  not  present 
in  measurable  amounts,  then  an  arbitrary  low  concentra¬ 
tion  must  be  assigned.  For  example,  if  the  major  con¬ 
stituents  are  present  at  the  10“-^  mol  L"'  range  and  Na 
was  not  detectable,  enter  an  arbitrary  low  concentration 
forNa  (e.g.,  10“**  mol  L“');  at  these  low  concentrations, 
Na  will  not  have  a  significant  impact  on  either  the  ionic 
strength  or  ion-pairing. 

If  there  are  constituents  not  explicitly  included  in  the 
program  that  are  present  in  significant  concentrations 
(e.g.,  NH4  or  NO3),  these  constituents  can  sometimes  be 
lumped  with  explicitly  recognized  constituents  as  input 
(e.g.,  NH4  with  K  [or  Na]  and  NO3  with  Cl).  Otherwise, 
the  program  will  have  to  be  modified  to  explicitly 
recognize  the  additional  constituents. 

The  program  was  written  to  analyze  laboratory  ex¬ 
periments  at  25°C.  As  a  consequence,  the  stability  con¬ 
stants  are  only  defined  for  25°C:  significant  departures 
from  this  temperature  will  require  program  modifica¬ 
tion.  For  example,  the  solubility  of  CO2  in  aqueous 
solution  is  about  15%  greater  at  20°Cthanat  25°C.  This 
solubility  difference  is  substantial  and  would  require 
program  modification  if  you  wanted  to  use  the  program 
at  2()''C. 

The  range  of  validity  of  the  lONPAlR  program  with 
respect  to  ionic  strength  (/)  is/ <0.1  M.  Monovalent  anil 
divalent  activity  coefficients  are  estimated  using  the 
Davies  equation  (eq  4),  which  is  considered  valid  for  all 
valency  ty  pes  up  to  /  =  0. 1  M,  with  an  error  not  exceeding 
2%  (Davies  1962).  Generally,  activity  coefficients  of 
neutral  molecules  deviate  from  1 .00  by  less  than  2-3% 
at/<0.1  M(Hamed  and  Owen  1958);  furthermore,  these 
deviations  may  be  both  greater  than  or  less  than  1.00. 
Given  the  small  magnitude  of  this  effect  on  neutrally- 
charged  constituents,  the  assumption  is  made  in  the 
lONPAlR  program  that  activity  equals  concentration 
( i.e.,  activity  coefficient  =  1 .00)  for  neutral  molecules. 

At  present,  the  computer  program  in  options  1  and  2 
assumes  that  measured  total  alkalinity  is  equivalent  to 
inorganic  carbon  alkalinity,  and  therefore  partitions  the 
alkalinity  amongst  the  inorganic  carbon  species.  If  this 
assumption  is  incorrect,  then  the  program  will  overesti¬ 


mate  the  true  inorganic  carbon  concentrations  and,  as  a 
consequence,  overestimate  the  true  ion  activity  products 
forCaC03.  However,  given  measurements  of  Pco,  tind 
pH  (option  3),  the  computer  program  will  calculate  the 
true  concentrations  of  inorganic  carbon  species,  which 
can  then  be  compared  to  the  measured  total  alkalinity. 
Recent  work  suggests  that  generally  total  alkalinity 
equals  inorganic  carbon  alkalinity  in  soils  (Suarez  1 990; 
Marion  et  al..  in  press  a);  but,  there  are  apparently 
exceptions  to  this  generality  when  total  alkalinity  does 
not  equal  inorganic  carbon  alkalinity  (Amrhein  and 
Suarez  1987,  Takkar  et  al.  1987,  Reddy  et  al.  1990). 


TESTING  CHEMICAL 
SPECIATION  MODELS 

The  validity  of  aqueous-phase  chemical  speciation 
models  can  be  examined  by  "overdetermining”  the  state 
of  these  aqueous  systems.  For  example,  given  any  two  of 
the  following  three  measurements — pH,  alkalinity  and 
Pco-,  — the  third  can  be  calculated.  For  lONPAlR's 
option  1  (pH  -f-  alkalinity),  Pco-,  can  be  calculated  from 
the  output  by 

Pco.  =  (H^)(HC03-)/K,  Kh 

=  (H^)(HC03-)/1.49x  10-**  (8) 

where  parentheses  refer  to  activities.  From  the  data  in 
Appendix  B  (option  1 ).  the  calculated  Pco,  (9-02  x  10" 
’  atm)  agrees  very  well  with  the  measured  PcOt  (8.96  x 
10 atm).  For  option  2  (alkalinity  and  Pco-|)>  calcu¬ 
lated  pH  is  given  by 

pH  =  -log  [(Pco,)  l  ‘^9x  10-**/(HC03-)|  .  (9) 

In  this  case,  the  calculated  pH  (7.29)  and  measured  pH 
(7.29)  agree  exactly  (Appendix  B).  And  last,  the  inor¬ 
ganic  carbon  alkalinity  can  be  estimated  from  option  3 
(pHaixI  Pco,)  (’V  summing  up  the  equivalent  concentra¬ 
tions  of  the  morganic  carbon  species.  In  this  case,  the 
calculated  inorganic  carbon  alkalinity  (3.35  x  10“-^moli; 
L  ')  is  in  agreement  with  the  measured  total  alkalinity 
(3.38  X  10"^  mole  )  (Appendix  B). 

In  this  particular  example,  the  calculated  and  mea¬ 
sured  quantities  (pH,  Pco,  ^  alkalinity)  were  in  excel¬ 
lent  agreement  as  were  tfie  charge  balances  (±  0.6  %, 
Appendix  B).  Over-determining  the  state  of  aqueous 
solutions  allows  one  to  check  the  internal  consistency  of 
experimental  measurements,  equilibrium  constants  and 
model  assumptions. 
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APPENDIX  A:  PROGRAM  LISTING 


1000  ! 

1010  !  Tins  PROGRAM  CALCULATES  SOIUriCN  ICNIC  OONCEMTRATICNS 
1020  !  WITH  OORBHCriCNS  FOR  ICN-PAIRS. 

1030  ! 

1040  OPTICM  NOLET 

1050  DIM  K(15),  X(20),  ACT  (20),  A(19,19) 

1060  DIM  B(19),  AC(20),  C(19),  AINV(19,19) 

1070  ! 

1080  !  SET  UP  SCREEM  AM)  PRINTER 
1090  ! 

1100  PAGE  =■  O.G 

1110  OPEN  #1:  SCREEN  0,1, 0,1 

1120  CPEN  #2:  PRINTER 

1130  PRINT  #2 

1140  ! 

1150  !  DEFINE  THE  STABILITY  OCNSTANTS 
1160  ! 

1170  DATA  7.43E-2,  7.04E-4,  4.90E-3,  8.50E-2,  5.75E-4 
1180  DATA  5.91E-3,  1.78,  5.39E-2,  1.41E-1,  1.78 
1190  DATA  5.39E-2,  2.01E-1,  4.68E-11,  4.33E-7,  3.43E-2 
1200  MAT  READ  K 
1210  ! 

1220  !  INEUT  DATA 
1230  ! 

1240  LINE  INPUT  PROMPT  "TITLE;  TITLES 

1250  PRINT  "CHOOSE  INPUT  CXTICN" 

1260  PRINT  "OPTION  1  =  ECD3T  +  PH" 

1270  PRINT  "OPTION  2  =  HOCI3T  +  002" 

1280  PRINT  "OPTICN  3  =  PH  +  002" 

1290  INPUT  PROMPT  "ENTER  NUMERIC  VALUE  OF  OPTICN  ":  OPT 
1300  INPUT  PROMPT  "ENTER  CA  OONC  (MOL/L)  ";  CALT 
1310  INPUT  PROMPT  "ENTER  tC  CCNC  (MOL/L)  ":  EGT 
1320  INPUT  PROMPT  "ENTER  K  OCMO  (MOL/L)  " :  KT 
1330  INPUT  PROMPT  "ENTER  NA  (XWC  (MOL/L)  ":  NAT 
1340  INPUT  PROMPT  "ENTER  S04  CCNO  (MOL/L)  ":  S04T 
1350  INPUT  PROMPT  "ENTER  CL  OONC  (MOL/L)  ";  CL 
1360  X{20)  =  CL 
1370  SEISCT  CASE  OPT 

1380  CASE  1 
1390  M>19 

1400  INPUT  PROMPT  "ENTER  PH  ":  PH 

1410  H-EXP(-2.3026*PH) 

1420  INPUT  PROMPT  "ENTER  ALKALINITY  (E)0/L)  ":  H003T 

1430  CASE  2 
1440  M>19 

1450  INPUT  PROMPT  "ENTER  AlEALINITY  (E)Q/L)  ":  HOOIT 

1460  INPUT  PROMPT  "ENTER  0)2  (ATM)  ":  C02 

1470  CASE  3 
1480  M>17 

1490  INPUT  PROMPT  "ENTER  PH  ":  PH 

1500  H=EXP(-2.3026*PH) 

1510  INPUT  PROMPT  "ENTER  002  (ATM)  ":  C02 

1520  EM)  SEliET 

1530  I4AT  REDIM  A(ND,M3),  B(ND),  C(M)),  AINV(MI,M)) 

1540  MAT  A  =  ZER 
1550  MAT  AINV  =  ZER 
1560  ! 

1570  !  ESTIMATE  INITIAL  IONIC  (XNCENTRATICNS 
1580  ! 

1590  SELECT  CASE  OPT 
1600  CASE  1 

1610  X(18)=HOO3T*0.95 

1620  X(19)=X(18)*K(13)/H 

1630  CTISE  2 

1640  X(18)'H003T*0.95 

1650  X(19)=X(18)*2*K(13)/(K(14)*K(15)*O02) 


S 


1660  CfiSE  3 

1670  X(18)-K(14)*K(15)*002/H 

1680  X(19)^(13)*X(18)/H 

1690  EM)  SEIfET 
1700  X(1)=0.85*CALT 
1710  X(2)-X(1)*X(18)/K(1) 

1720  X(3)-X(1)*X(19)/K(2) 

1730  X(4)=X(l)»SOST/K(3) 

1740  X(5)-0.85*M3T 
1750  X(6)=X(5)*X(18)/K(4) 

1760  X(7)=X(5)*X(19)/K(5) 

1770  X(8)=X(5)*S04T/K(6) 

1780  X(9)=0.97*KT 

1790  X(10)=X(9)*X(18)/K(7) 

1800  X(11)=X(9)  *X(19)/K(8) 

1810  X(12)=X(9)*X(17)/K(9) 

1820  X(13)=0.97*NAT 

1830  X(14)=X(13)‘X(18)/K(10) 

1840  X(15)=X(13)*X(19)/K(11) 

1850  X(16)=X(13)*X(17)/K(12) 

1860  X(17)=SO4T*0.85 
1870  ! 

1880  !  ESTIMATE  THE  IONIC  STRENGTH  AM)  THE  MONO-  AM)  DI-VALEMT 
1890  ACTIVITY  OOEFFICIEMrS  USINS  THE  DAVIES  EQUATICN. 

1900  ! 

1910  ICNSTFK).5*(X(2)+X(6)+X(9)+X{11)+X(12)+X(13)+X(15)+X(16)+X(18)+X(20)) 
1920  laiSTT^ICNSTR  +  2.0*  (X(l) +X(5) +X(17) +X(19) ) 

1930  FACroR-SQR(ICNSTO)/(1.0+SC2R(IONSTR))-O.3*ICNSTR 
1940  ACl=EXP(-1.1720*F«nOR) 

1950  AC2=EXP(-4.6881*FACrOR) 

1960  ! 

1970  I  DEFINE  THE  EXEMNTS  OF  TOE  PARTIAL  DERIVATIVE  (A)  AM)  THE 
1980  !  DIFEXREM3;  (B)  MATRICES  FOR  THE  NEWICN-RAPHSCN  AIGORITIW 
1990  ! 

2000  A(l,l)=-1.0 
2010  A(l,2)=-1,0 
2020  A(l,3)=-1.0 
2030  A(l,4)=-1.0 
2040  A(2,5)=-1.0 
2050  A(2,6)=-1.0 
2060  A(2,7)=-1.0 
2070  A(2,8)=-1.0 
2080  A(3,9)=-1.0 
2090  A(3,10)=-1.0 
2100  A(3,ll)=-1.0 
2110  A(3,12)=-1.0 
2120  A(4, 13)=-1.0 
2130  A(4, 14)=-1.0 
2140  A(4,15)=-1.0 
2150  A(4,16)=-1.0 
2160  A(5, 17)=-1.0 
2170  A(5,4)=-1.0 
2180  A(5,8)=-1.0 
2190  A(5,12)=-1.0 
2200  A(5,16)=-1.0 
2210  AK1=K(1)/AC2 
2220  AK2=K(2)/AC2'2 
2230  AK3=K(3)/AC2'2 
2240  AK4=K(4)/AC2 
2250  AK5=K(5)/AC2"2 
2260  AK6=K(6)/AC2''2 
2270  AK7=K(7)/A£;i''2 
2280  AX8=K(8)/«2 
2290  AK9=K(9)/AC2 
2300  AK10=K(10)/AC1"2 
2310  AKll-K(ll)/«2 
2320  AK12-K(12)/AC2 
2330  A(6,l)-X(18) 

2340  A{6,2)— AKl 
2350  A(7,1)=X(19) 
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2360  A(7,3)^«<2 
2370  A(8,1)-X(17) 

2380  A(8,17)“X(1) 

2390  A(8,4)— AK3 
2400  A(9,5)“X(18) 

2410  A(9,6)=-AK4 
2420  A(10,5)=X(19) 

2430  A(10,7)»'AK5 
2440  A(11,5)“X(17) 

2450  A(11,17)-=X{5) 

2460  A(11,8)='AK6 
2470  A(12,9)“X(18) 

2480  A(12,10)=-AK7 
2490  A(13,9)=X(19) 

2500  A(13,11)=-AK8 
2510  A(14,9)=X(17) 

2520  A(14,17)=X(9) 

2530  A(14,12)=-AK9 
2540  A(15,13)=X(18) 

2550  A(15, 14)=-AK10 
2560  A(16,13)=X(19) 

2570  A(16,15)=-AK11 
2580  A(17,13)=X(n) 

2590  Aa7,17)=X(13) 

2600  A(17,16)=-AK12 
2610  SEIiCT  CASE  OPT 
2620  CASE  1  TO  2 
2630  A(18,18)=-1.0 

2640  A(18,2)=-1.0 

2650  A(18,6)=-1.0 

2660  A (18, 10) =-1.0 

2670  A(18,14)=-1.0 

2680  A(18,19)=-2.0 

2690  A(18,3)— 2.0 

2700  A(18,7)— 2.0 

2710  A(18,ll)=-2.0 

2720  A(18,15)=-2.0 

2730  A(6,18)=X(1) 

2740  A(7,19)-X(l) 

2750  A{9,18)=X(5) 

2760  Ado,  19)  =X  (5) 

2770  Ad2,18)=X(9) 

2780  A(13,19)=X(9) 

2790  A  (15, 18)  =X  (13) 

2800  A(16,19)-X(13) 

2810  B(18)=H303T-X(18)-2.0*X(19)-X(2)-2.0‘X(3)-X(6)-2.0*X(7)-X(10) 

2820  B(18)=B(18)-2.0*X(11)-X(14)-2.0*X(15) 

2830  SEIBCT  CASE  OPT 

2840  CASE  1 

2850  AK13=K(13)*AC1/(AC2*H) 

2860  A(19,19)=1.0 

2870  A(19, 18)=-AK13 

2880  B(19)=X(19)-AK13‘X(18) 

2890  case  2 

2900  AK14=K(14)  ‘KdS)  *002*^72/ (K(13)  •AC1''2) 

2910  A(19,18)=2.0‘X(18) 

2920  A(19, 19)=-AK14 

2930  B(19)=X(18)'2-AK14*X(19) 

2940  Q®  SEIECT 

2950  CASE  3 

2960  X(18)=K(14)*O02*K(15)/(H*ACl) 

2970  X(19)-K(13)*X(18)*AC1/(H*AC2) 

2980  EM)  SElfXTT 

2990  B(1)=CALT-X(1)-X(2)-X(3)-X(4) 

3000  B(2)>+CT-X(5)-X(6)-X(7)-X(8) 

3010  B(3)-KT-X(9)-X(10)-X(11)-X(12) 

3020  B(4)-NAr-X(13) -X(14)  -X(15)  -X(16) 

3030  B(5) -S04T-X(17) -X(4) -X(8) -X(12) -X(16) 

3040  B(6)-X(1)*X(18)-AK1*X(2) 

3050  B(7)-X(1)*X(19)-AK2*X(3) 
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3060  B(8)=X(1)*X(17)-AK3»X(4) 

3070  B(9)“X(5)*X(18)-AK4*X(6) 

3080  B(1C)=X(5)»X(19)-AK5*X(7) 

3090  B(n)«X(5)*X(17)-AK6*X(8) 

3100  B(12)-X(9)*X(18)-AK7*X(10) 

3110  B(13)-X(9)*X(19)-AK8*X(11) 

3120  B(14)“X(9)*X(17)-AK9*X(12) 

3130  B(15)=X(13)  *xa8)-AK10*X(14) 

3140  B(16)=X(13)*X(19)-AK11*X(15) 

3150  B(17)=X(17)*X(13)-Aia2*X(16) 

3160  FOR  1=1  TO  ID 
3170  B(I)— 1.0*B(I) 

3180  lEXT  I 
3190  ! 

3200  !  SOLVE  THE  SET  OF  LINEAR  EQUATIONS  '0=B/A) 
3210  ! 

3220  MAT  AW  -  INV(A) 

3230  MAT  C  =  AINV*B 
3240  ! 

3250  !  QECK  FOR  SOLUTICN  OCNVERGINCE 
3260  ' 

3270  FOR  1=1  TO  tD 
3280  X(I)=X(I)-K;(I) 


3290 

3300 

3310 

3320 

3330 

3340 

3350 

3360 

3370 

3380 

3390 

3400 

3410 

3420 

3430 

3440 

3450 

3460 

3470 

3480 

3490 

3500 

3510 

3520 

3530 

3540 

3550 

3560 

3570 

3580 

3590 

3600 

3610 

3620 

3630 

3640 

3650 


NEXT  I 

for  1=1  TO  ND 

PCEN=(C(I)/X(I))*100.0 
PCE»=ABS  (PON) 

IF  PCIK>0.1  THEN  GO  TO  1910 
NEXT  I 

!  PRINT  RESULTS 


PRINT  #2 
PRINT  #2 

PRINT  #2:  TAB(7);  "INPOT" 
PRINT  #2;  TAB (8);  TITLES 
PRINT  #2,  USQC  " 

SELECT  CASE  OPT 
CASE  1 

PRINT  #2,  USING  " 
CASE  2 

PRINT  #2,  USING  " 
CASE  3 

PRINT  #2,  USING  " 


CAT=#. ##""""  KT=#.#r"'''  NAT=l.#r CALT,M3T,iCT,NAT 


iX»T=#.  CL=#.##^''^  ALK=#.##'' 

S04T=# .  CL=#.##'^''''  AIEC=#.##" 


S04'IH1.## -  (3/=#.## -  PH=*.###  C02=#.##-"' 

H003I^X  (2) +2*X  (3) +X  (6) +2*X  (7) +X  (10) +2»X  (ID +X  (14) +2*X  ( 16) 

K»3T=)CD3T+X(18) +2 .0*X(19) 

EID  SEIECT 

C7)TICN=2 . 0‘  (CALT-rtCT)  +KT+NAT 
ANICN=2 . 0*S04TKX+H003T 
PEK3Nr=(  (CATICN-ANICN)/CATICN)  *100.0 
PRINT  #2 

PRINT  #2;  TAB  (7);  "COTPOT" 

PRINT  #2,  USING  "  TOTAL  CATIONS  TOTAL  ANICNS=#.## - :  CATION, ANICN; 

PRINT  #2,  USING  "  PERCENT  DIFF=  +♦»#.##”:  PERCENT 

PCa)'-K2 

AC(2)-AC1 

AC(3)-I.0 

AC(4)=1.0 

fC(.5)=N2 

«:(6)=PC1 


PH=#.###”:  S04T,a,HC03T,PH 

S04T,CL,H003T,C02 
S04T,CL,Pri,C02 


3560  a; (7) =1.0 
3670  AC(8)=1.0 
3680  AC(9)=A.'1 
3690  AC (10) =1.0 
3700  AC(ni=ACl 
3710  AC(12)=A;1 
3720  AC(13)=AC1 
37  30  AC(14)=1.0 
3740  AC(15)=AC1 
3750  a:(16)=AC1 
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3760  AC(17)-=flC2 

3770  PCim-fCl 

3780  AC(19)-AC2 

3790  «:(20)=«:i 

3800  FOR  1=1  TO  20 

3810  lCr(l)=X(l)*AC(I) 

3820  NEXT  1 


3830 

PRINT  #2:  " 

SPECIES  GONG.  ACT. 

SPECIES  oac.  ACT. 

" 

3840 

PRINT  #2,  USING  " 

CA=  # 

CAH003=#.##^ - 

I.##'"-''" 

X(1),ACT(1),X(2),ACT<2) 

3850 

PRINT  #2,  USING  " 

CACD3=#.  ♦ 

CAS04= 

X(3),ACT(3),X(4),ACT(4) 

3860 

PRINT  #2,  USING  " 

M>  ♦.I#''""''  t 

MGHOOSH.##"^"" 

t.tt' - 

X(5),ACT(5),X(6),ACT(6) 

3870 

PRINT  #2,  USING  " 

M3003=-#.##'"^"  « 

AA  A/VAA 

MGS04=  f.##^ - 

I.##"""'" 

X(7),ACT(7),X(8),ACT(8) 

3880 

PRINT  #2,  USING  " 

K=  * 

KHC03=  #.#» - " 

X(9),ACT(9),X(10),ACT(10) 

3890 

PRINT  #2,  USING  " 

KC03=  « 

AAAAA A 

KS04= 

X(ll)  ,ACT(11),X(12),ACT(12) 

3900 

PRINT  #2,  USING  " 

NA=  ♦.##-''""  # 

NAHOO>#.## - 

#.#1"^ - - 

X(13),ACT(13),X(14),ACT(14) 

3910 

PRINT  #2,  USING  " 

NAO03H».##''"""  # 

NAS04- 

X(15),ACT(15),X(16),ACT(16) 

3920 

PRINT  #2,  USING  “ 

S04*  « 

H003=  I.##""''" 

#.##"^ - - 

X(17),ACT(17),X(18),ACT(18) 

3930 

PRINT  #2,  USING  " 

003=  # 

CL=  #.♦# - 

X(19),ACT(19),X(20),ACT(20) 

3940  CAS=a7r(l)*ACr(17) 

3950  CAS=-IO310(CAS) 

3960  A2TPRCD=*Cr(l)*ACT(19) 

3970  P1AP=-113G10  (flCTPBCD) 

3980  PRINT  #2,  USING  "  IONIC  STRENCTWI.##^^-'''  PIAP(CB003)  =#».###";  ICNSTR, PIAP; 

3990  PRINT  #2,  USING  "  PIAP(CAS04)=##.###":  CAS 

4000  PAGE=Pft3;+1.0 

4010  IF  PAGE=2.0  THIN 

4020  PRINT  #2:  CHRS(12) 

4030  PAOX.O 
4040  ELSE 
4050  PRINT  #2 
4060  DC  IF 
4070  ! 

4080  !  CHECK  FOR  MORE  DATA 
4090  ! 

4100  INPUT  PRCMPT  "IS  THERE  MORE  DATA  (Y(l)  OR  N(0))?":  MORE 
4110  IF  MORE  =  1  THEN  GO  TO  1240 
4120  IF  PAOX)  THEN  GO  TO  4140 
4130  PRINT  #2;  CHR$(12) 

4140  END 
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APPENDIX  B:  SAMPLE  INPUT  AND  OUTPUT 


INPUT 

ALTAMa^  -  41  DAYS  -  OPTICN  1 

CAT=1.07e-02  M3T=3. 52^03  KT=1.13e-04  NAT=6.09e-04 
SO4T=1.28e-02  CL=5.08e-05  ALK=3.38e-03  PH=7.290 

OUTPUT 

TOTAL  CATIONS=2.92e-02  TCOAL  ANIGNS=2.90e-02  PEBCEOT  DIFF=  +  .45 


SPECIES  OC*4C. 

ACT. 

SPECIES  0340. 

ACT. 

CA=  7.41e-03 

3.60e-03 

CAHOO3=1.52e-04 

1.27e-04 

CAOO>1.22e-05 

1.22e-05 

CAS04=  3.13e-03 

3.13e-03 

M3=  2 . 57e-03 

1.25e-03 

MC31003=4.61e-05 

3.85e-05 

MGOO3=5.19e-06 

5.19e-06 

M3S04=  8.99e-04 

8.99e-04 

K=  l.lOe-04 

9.15e-05 

KH003=  1.35e-07 

1.35e-07 

K003=  4.86e-09 

4.05e-09 

KS04=  3.31e-06 

2.76e-06 

NA=  5.96e-04 

4 . 97e-04 

NAH003=7. 31^07 

7.31e-07 

NA003=2 . 64e-08 

2.20e-08 

NAS04=  1.26e-05 

1.05e-05 

S04=  8.76e-03 

4.26e-03 

HQ03=  3.14e-03 

2.62e-03 

003=  4 . 92e-06 

2.39e-06 

CL=  5.08e-05 

4.24e-05 

KMC  STPENGrH=3.95e-02  PIAP(CA003)=  8.065  PIAP(CAS04)=  4.815 


INPUT 

ALTAMM  -  41  DAYS  -  (MION  2 
CAT=1.07e-02  M3T=3.52e-03  K]>=1.13e-04  NAT=6.09e-04 
SO4T=1.28e-02  CL=5.08e-05  AIi<=3 . 38e-03  CO2=8.96e-03 

OUTPUT 

TOTAL  CATia4S=2 . 92e-02  TOTAL  ANICNS=2.90e-02  PERCENT  DIFF=  +  .45 


SPECIES  OONC. 

ACT. 

SPECIES  OCNC. 

ACT. 

CA=  7.41e“03 

3.60e-03 

CAfKX)3=1.52e-04 

1.27e-04 

CA003=1.23e-05 

1.23e-05 

CAS04=  3.13e-03 

3.13e-03 

MG=  2.57e-03 

1.25e-03 

M3i003=4.61e-05 

3.85e-05 

MXO3=5.23e-06 

5.23e-06 

M3S04=  8.99e-04 

8.99e-04 

K=  1 . lOe-04 

9.15e-05 

KH003=  1.35e-07 

1 . 35e-07 

KOO>  4.90e-09 

4.09e-09 

KS04=  3.31e-06 

2.76e-06 

NA=  5.96e-04 

4 . 97e-04 

NAH003=7.31e-07 

7.31e-07 

NA003=2 . 66e-08 

2.22e-08 

NAS04=  1.26^05 

1.05e-05 

S04=  8.76e-03 

4.26e-03 

H003=  3.14e-03 

2.62e-03 

003=  4 . 96e-06 

2.41e-06 

CL=  5 . 08«=-05 

4 . 24e-05 

IONIC  STRENGTH=3 . 95e-02  PIAP(CA003)=  8.062  PIAP(CAS04)=  4.815 


li 


INPUT 

ALTAMCNT  -  41  DAYS  -  CPTION  3 

CAT=1.07&-02  MGT=3.52e-03  KT=1.13e-04  NAT=6.09e-04 

SO4T=1.28e-02  CL=5.08e-05  PH=7.290  002=8. 96e-03 

OUTPUT 

total  CATiaiS=2.92e-02  TOTAL  ANIONS=2.90e- 02  PERCENTT  DItT=  + 


SPECIES  COOC. 

ACT. 

SPECIES  OONC. 

ACT. 

CA.=  7.41e-03 

3.60e-03 

CAHOO3=1.51e-04 

1.26e-04 

CAOO3=1.21e-05 

1.21e-05 

CAS04=  3.13^03 

3.13e-03 

M3=  2.57e-03 

1.25e-03 

MC51OO3=4.57e-05 

3.81e-05 

M3Cr)3=5.14e-06 

5.14e-06 

MSS04=  8.99e-04 

8.99e-04 

K=  1 . lOe-04 

9.15e-05 

KH003=  1.33e-07 

1.33e-07 

KO03=  4.81e-09 

4 . 02e-09 

KS04=  3.31e-06 

2.76e-06 

NA=  5 . 96e-04 

4 . 97e-04 

NAH003=7.25e-07 

7.25e-07 

NA003=2 . 62e-08 

2.18e-08 

NAS04=  1.26e-05 

1.05e-05 

S04=  8.76e-03 

4.26e-03 

H003=  3.11e-03 

2.59e-03 

003=  4.87e-06 

2.37e-06 

CI^  5.08e-05 

4 .24e-05 

IONIC  STRENGTH=3.95e-02  PIAP(CA003)=  8.069  PIAP(CAS04)=  4.815 
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REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMBNo.  0704-0188 


PuDHC  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existtrig  data  sources  gathering  and 
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